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Abstract 

In this paper we consider a new nonUnear mathematical model for 
corneal topography formulated as two-point boudary value problem. 
We derive it from first physical principles and provide some mathe- 

Smatical analysis. The existence and uniqeness theorems are proved 
, as well as various estimates on exact solution. At the end we fit the 

simplified model based on Modified Bessel Function of the First Kind 
with the real corneal data consisting of matrix of 123 x 123 points and 
obtain an error of order of 1%. 

1 Introduction 

^ Sight is the most important sense that we posses. It enables us to truly 

perceive the outer world. It is crucial to understand the mechanics of vision 
K*" in order to treat various diseases that may occur and disturb regular seeing. 

^ The eye's main part responsible for about two-thirds of refractive power is 

^ the cornea. It is one of the most sensitive parts of human body and its various 

irregularities can cause many seeing disorders. Common diseases like myopia, 
hyperopia or astigmatism have their origin in some abnormalities in corneal 
geometry. Thus, precise knowledge of corneal shape is very important. 

We can distinguish several mathematical models of cornea which are 
presently in use. Most commonly, cross-section of cornea is described by 
some conic section [T] (for example, parabola or ellipse). The drawback of 
these models is that they have almost no physical motivation or derivation 
behind them. Probably the most appropriate model of corneal geometry 
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would be done with the use of shell theory [2j. These models usually are 
quite complex and hence difficult to analyse. There are also models that use 
Zernike polynomials These are certain orthogonal polynomials that are 
used to model abberations in cornea or lens. 

In this paper we present another model of corneal geometry based on a 
nonlinear membrane equation. Assuming radial symmetry we reduce it to a 
one- dimensional case and then provide some estimates on the exact solution. 
When fitting with data we use its simplified form and find that mean error 
is of order of a few per cent. 

2 Derivation 

Cornea is a biological structure which is similar to a shell. Its diameter is 
about 11mm and its thickness is approximately 0.5mm in the center and 
0.7mm in the peripheral part jl]. Thanks to its structural parameters cornea 
provides a great protection for more sensitive, inner parts of the eye. We will 
distinguish only interior and exterior surface of cornea although it is made of 
five layers. In order from the most exterior they are as follows: epithelium, 
Bowman's layer, stroma, Descement's membrane and endothelium. Stroma 
is the thickest layer sometimes identified with cornea itself. 

We would like to derive an equation describing the steady state shape of 
cornea's surface (either interior or exterior). We will describe this geometry 
by a function h = h{x,y). Assume that cornea is a two-dimensional thin 
membrane, that is, it has no bending or twisting moments and its tension 
T is constant at every point. Moreover, impose a pressure P acting from 
below and normally to the surface. To model cornea's elastic properties we 
propose an additional restoring force which is proportional to the defiection 
h. Finally, let Q be the domain in xy plane on which this membrane is 
situated. 

It can be shown (see [3]) that a steady state membrane equation is of the 
form 

- TAh = F, (1) 

where F is the net force density (pressure) acting on the membrane. In our 
case the equation will be 

-TAh + kh= ^ =. (2) 

Vi + liv/^ir 

This is a second order nonlinear partial differential equation. The square 
root on the right-hand-side comes from the projection of normal vector onto 
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the vertical z-axis. For the boundary conditions we impose 



h = on dQ. (3) 

Appearing constants have following dimensions: [T] = — = MT~'^, [k] = 
4 = ML-'^T-^ and \P] = ^ = ML~^T'^, where M, L^T are dimensions 
of mass, length and time respectively. Also, let i? be a typical linear size of 
the cornea (for example its radius). Taking 

h* = - X* = - (A) 
we transform ^ into a nondimensional equation (dropping asterisks for clar- 

ity) 

— A/i + ah = — = on f2, (5) 

^Jl + \\Wh\\' 

where we denoted a := and b := From physical reasoning we allow 
only positive values of a and b. We will see in Section |4] that these parameters 
are both of order of unity. 



3 Model Analysis 

From now on Q will be an unit circle in the xy-plane. We assume that 
cornea has radial symmetry, thus h depends only on radius r and in polar 
coordinates h = h{r). Rewriting ([s]) we get 

I d f dh\ ^ b 

[r—\ +ah= , = 0<r<l 6 

rdr\ dr ) ^TTTi^ ~ ~ 

with following boundary conditions 

/i'(0) = 0, h{\) = 0. (7) 

The condition h'{0) = guarantees absence of singularities and smoothness 
at the origin. 

We transform ^ into an integral equation 

h{r) = (^,{r) j\v,{t)P{h\t))dt + v^{r) j\vo{t)P{h'{t))dt^ , 

(8) 
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where Vq, Vi and P are defined by 

Vo{r) := Io{y/ar), Vi{r) := Io{^/a)Ko{^/ar) - lQ{y/ar)KQ{y/a), 

1 (9) 



Fix) 



and J,^, Ki, {v G M) are i/-order modified Bessel functions of first and second 
kind respectively. From elementary properties of Bessel functions (see [6]) 
we can deduce the following estimates 

t;o(r)>0, t;i(r)>0, v'^{r) = ^h{^r) > 

v[{r) = -v^ (/o(v^)iri(v^r) + h{,/Er)Ko{y^)) < 0. 

Thus, we have h > 0. Also, it can be shown that 

i- {-rv[{r)) = -ar {h{^)K^{^r) - K,{^)h{^r)) < 0, (11) 

hence the function —rv'{r) is nonincreasing. Since Ky[z) ~ }^V{u){j^z)^'^ 
we have \mir^Q+ {—rv'{r)) = iQ^y/a). Moreover, notice the following useful 
formulas [B]: 

Uz)K,+^{z) + h+i{z)K,{z) = - (12) 



and 



zK,{z) = {zK,{z)) , zioiz) = ^ {zh{z)) , 

-f/o(z)=/i(z), iro(^) = -i^i(^). 
02; az 



(13) 



Finally, it is easy to calculate that 



b 



h'{r) = —^iv',{r) tvi{t)P{h'{t))dt + v[{r) tv,{t)P{h' {t))dt 

(14) 

The question about existence and uniqueness of solution to ^ is answered 
by the following theorem which uses Picard-Lindeloeff iterations. 

Theorem 1. Define the approximation sequence by 

a \ /o(va) / 

hnir) := (^voir) tvi{t)P{h'^_^{t))dt + Vi{r) tvQ{t)P{h'^_^{t))dt 

(16) 
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where n = 1,2, ... and < r < 1. The sequence hn is uniformly convergent 
to the unique solution of boundary value problem @-([^ provided that 

, 3v^ v/^/o(v/a) 

2 /i(^)(2/oV^-l)- ^ ^ 



Proof. First, observe that P is Lipschitz continuous: \P{x) — P{y)\ < M \x — y\, 
where M = We can write (8) and (14) as h{r) = F{r,t)P{h'{t))dt 

and h'{r) = G{r,t)P{h'{t))dt, 



wnere 



F{r,t) 



{tvo{r)vi{t)x[i,r]i't) + tvo{t)vi{r)x[o,r]it)) 



G{,r,t) := {tv'Q{r)vi{t)x[i,r]{t) + tvo{t)v[{r)x[o,r]{t)) . 

It is an easy calculation using mJ-(Q3f to show that /J \F{r, t)\dt<Q and 
jl \G{r,t)\dt < R, where 



R ■-- 



(19) 



(2/ov^-l). 



^/aIo{^/a) 

Let hn{r) = Hq^t) + J27=i {hi{r) - K^iir)) and similarly h'^{r) = h'^ir) + 
X^ILi {Ki^) have to show that these series are uniformly 

convergent when n — )■ oo. By estimate 

l^'i(r) - K{r)\ < f \G{r,t)\ \Pm)) - P(0)| 



\G{r,t)\ \h'Q{t)\dt < MR^, 



< M 



and generally 

Kir) - K,_,{r)\ < f \G{r,t)\ \P{h'^_,{t)) - P{h'^_,{t))\dt 
Jo 

<M [ \Gir,t)\K^,it)-K^^,{t)\dt<MR\K_,-h:_,l 
Jo 



we get \K - < MR \\K-i - K 

supremum norm. Moreover, 



^ < R{MRY, where ||.||^ is a 

\K{r) - /i„_i(r)| <M f \F{r,t)\ %_,{t) - h'^_,{t)\dt < Q^MRY, 

Jo 



that is \\hn — hn-i\\^ < Q{MR)^. By the assumption MR < 1, so the 
dominant number series are uniformly convergent. From uniform Cauchy- 
sequence argument, hn ^ f and h'^ — /' for some /. The equahty f = h 
follows from continuity of P and (16). This ends the proof. □ 



Imposing more restrictions on b we can prove some estimates for the 
solution to (IgI). 



Lemma 2. Assume that b < , - ^^,°^/^\ ^ , then 



1 A 



h' <h' <0. 



(20) 



Proof. Take b as in the assumptions. By (10)-(14) and since P < 1 we have 
the estimate 

"1 



(21) 



h\r) < J^^<ir) tv,{t)P{h'{t))dt 
bhiJar 



< 



Now, by (11) and our assumption on b we deduce that 

bh{^a) 



h'(r) < 



2/o(v/S) - 1 



( r\ (^o(v^) - 1) < ^ . . 
Next, by similar reasoning we have 



J, r-Mr) f tvo{r)P{h\t))dt < - ^!^lyj-. rv[{r) 
hWa) Jo y/aIo{y/a) 

^ bh{^) ^ v/2/o(v/^) - 1 

- /o(v^)-l 



We have shown that \h'{r)\ < ^f^l^, thus by (|9|) 



C < P{h'{t)) < 1, 



(22) 



(23) 



(24) 



where C = 1 — (0 < C < 1). Finally, from (11) and (24) we deduce 
that 



h'ir) < 



v'oir) / tvi{t)dt + Cv[{r) / tvo{t)dt 

-'OlVO-j V Jr Jo 



^Jlip^ ("(1 _ C)r (/o(v^)Ki(v^r) + Ko{V^)h{V^r)) - < 0. 

(25) 



Similarly, estimating from below we can prove 



-h'{x) < - !^ fgfo(r) / tvi{t)dt + v[{r) [ tvo{t)dt 
-'olVaj V J J- Jo 



Io{y/a) V Va 



^bhiV^(^_ 1 



■v/a/o(i/a) 



Kir). 



(26) 

□ 



Therefore we have proved our claims. 

Corollary 1. Let b < ^ ■ Solution of the problem (6)-(l) is a 

positive, nonincreasing function h with 



Ahi <h<hQ 



where ho and hi are defined by (15)-(16), while 

A= — 



1+2- 



h'il] 



(27) 



(28) 



Proof. Lemma |2] gives us the fact that h is nonincreasing. Moreover, we have 
P ^^2 — jp(y^) j ^o) — — 1- Simple calculus allows us to deduce that 

AP (h^) < P ^^2 — /q(^) j h'o^. Now using the equation for h (8) we can 
conclude the proof. □ 



Various estimates on h are depicted on Figure [Tj The convergence of 
hn is very fast. In fact, only after 4 iterations we have {{h^ — h^W^ ^ 10~^. 
Although Picard's iterations can be very demanding on computer power they 
converge very rapidly hence the number of iterations can be kept small. 

To provide validity of Theorem [T] and Corollary [T] we have to restrict the 
values of b. Nevertheless, we will see in the next section that values of (a, b) 
calculated from the data lie in the admissible region. Figure [2] shows the 
assumed restrictions on b. 



4 Data Fitting 

In this section we present the results of data fitting. Instead of exact solution 
of ([6| we will use its zeroth-order approximation ho which is defined by (15). 
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Figure 1: Estimates on /i for a = 6 = 2. Dashed line: ho (top) and Ahi 
(bottom). Solid line: h^. 




Figure 2: Values of (a, 6) assumed in T heorem [I] and Corollary [T] Solid line: 



^ Dashed line: 

2 /iv^(2/ov^-l) 



^ ^/2Jo(v/^)-l 



The use of ho rather than h is equivalent to making the assumption that the 
pressure acts in the z-direction as apposed to the normal direction. Moreover, 
ho has a simple analytic form and hence is much easier to use in applications. 

We will determine constants a and b by taking two measurements: cornea's 
maximum deflection and central radius of curvature. These two values are 
readily measurable quantities and allow us to make a fit on the whole domain. 
Suppose we know maximum deflection ho{0) and central radius of curvature 

p(0), where p{x) = . Because /i'o(O) = we have p(0) = 

thus 
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where in calculating the second derivative we used 21[{z) = Io{z) + hiz) 
(see j6]). Finally, using (15) we deduce that a is the solution of the following 
equation 



1 



/io(0)p(0)a - /o(v^) + 1 = 0. 



(30) 



For two datasets consisting of 123 x 123 points representing interior and ex- 
terior surfaces of the cornea we calculated that ~ 2.07883, 6j„ ~ 2.76741, 
ttex ~ 1.94398 and bex ~ 2.27534. It is evident that these points are within 
the admissible region depicted on Fig. |2} 

To make the fitting more accurate we can resign from assumed radial 
symmetry by making a slight perturbation in the domain Q. Specificly, as 
Q we will use ellipse instead of a circle. The eccentricity of this ellipse can 
be calculated from the data, for example, by measuring the eccentricity of 
level curves of the corneal shape. In order to reshape our Q into an ellipse 

we have to change r = a/x^ + into r = where Ri and R2 are 

semiaxes of the ellipse in the domain. 

Mean fitting errors for interior and exterior surfaces are respectively 
0.035mm (3.6%) and 0.014mm (1.4%) and squared eccentricities are —0.0214 
and 0.0234. The interior surface is thus slightly prolate while exterior is 
oblate. It is evident that deviations from circular shape are very minute. 
The magnitude of absolute fitting error is shown on Fig. [3j 



Absolute Error for fitting (Interior Surface] 





Figure 3: Absolute fitting errors for interior (left) and exterior (right) sur- 
faces. 

Another important feature of cornea's geometry is its curvature. For 
measuring radius of curvature optometrists use the following quantity known 
as axial distance ([7]) 
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In the case of surface of revolution, normal vector lies in the meridional plane 
so d is the distance along the normal vector from the surface point to the 
axis of revolution. The absolute error in fitting d calculated for our is 
depicted on Fig. |4| Mean errors are approximately 0.146mm (2.1%) and 
0.14mm (1.7%) for interior and exterior surfaces respecively. 



Absolute error for axial curvature (Interior Surface) Inim] Absolute error for axial curvature (Exterior Surface] [mm] 




Figure 4: Absolute errors in axial distance for interior (left) and exterior 
(right) surfaces. 



5 Conclusion and discussion 

We have proposed a nonlinear mathematical model of corneal geometry based 
on a thin membrane equation. The existence and uniqueness of the solution 
was proved provided that one of the parameters was bounded in respect the 
the other. Various estimates (Lemma [2] and Corollary [T]) showed that the 
solution behaves according to common physical intuition. Approximating se- 



quence (16) converges very rapidly to the exact solution what can be deduced 
from its definition and numerical experiments. 

We have fitted zeroth-order approximation to the two meshes of 123 x 
123 points representing interior and exterior surfaces of cornea. We have 
slightly perturbed its radial symmetry according to the measurable quanti- 
ties. The mean fitting error was of order of a few per cent which provided 
a good fit. The error was slightly smaller for exterior than interior surface. 
Also, the function known as axial distance representing cornea's radius of 
curvature was in a good accordance with the data. Again, the magnitude of 
errors were of order of a few per cent what allows us to conclude that this is 
a good result for a model of that simplicity. 
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